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Abstract 

A new calculation of the atmospheric fluxes of cosmic-ray hadrons and muons in 
the energy range 10-10 5 GeV has been performed for the set of hadron production 
models, EPOS 1.6, QGSJET 11-03, SIBYLL 2.1, and others that are of interest 
to cosmic ray physicists. The fluxes of secondary cosmic rays at several levels in 
the atmosphere are computed using directly data of the ATIC-2, GAMMA exper- 
iments, and the model proposed recently by Zatsepin and Sokolskaya as well as 
the parameterization of the primary cosmic ray spectrum by Gaisser and Honda. 
The calculated energy spectra of the hadrons and muon flux as a function of zenith 
angle are compared with measurements as well as other calculations. The effect of 
uncertainties both in the primary cosmic ray flux and hadronic model predictions 
on the spectra of atmospheric hadrons and muons is considered. 
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1 Introduction 



A comparison of the calculated spectra and zenith-angle distributions of atmo- 
spheric hadrons and muons with the results of measurements makes it possible 
to solve related problems: (i) indirect research of the high-energy behavior of 
the primary cosmic ray spectrum and composition provided that cross sections 
of hadron-nuclear interactions are known, and (ii) a study of the high-energy 
hadron-nucleus interactions provided that energy spectra and composition of 
primary cosmic rays or spectra of secondaries are measured with a satisfactory 
accuracy. 

In this work we analyze the possibility to discriminate indirectly hadronic in- 
teraction models leaning on recent measurements of hadron and muon fluxes. 
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We present new calculations of atmospheric cosmic-ray fluxes in the range 10- 
10 5 GeV performed with current models extensively used in the simulation of 
high-energy cosmic ray propagation through the atmosphere, QGSJET01 [lj, 
QGSJET II [2,3J, SIBYLL @], NEXUS 0, EPOS [6J and Kimel and Mokhov 
(KM) |7j. Numerical results are obtained using directly data for the primary 
cosmic ray spectra and composition of the ATIC-2 [8] and GAMMA [9] ex- 
periments. Alternatively we use recently proposed model by Zatsepin and 
Sokolskaya (ZS) [10] as well as the parameterization by Gaisser and Honda 
(GH) [TT]. The comparison of the results of these calculations with recent di- 
rect measurements of high-energy atmospheric muons [1"2"1|13II14|[T5] and with 
predictions of other authors allow to evaluate uncertainties caused by statisti- 
cal errors of the primary spectrum data as well as to discriminate them from 
those originated from hadronic models. 

The calculations are based on the method originally developed for solving 
problems of neutrino transport in matter [16J and modified subsequently [T7] 
to apply for the transport problem of the nucleon component of cosmic rays 
in the Earth's atmosphere. The solution of the nucleon-meson cascade equa- 
tions was stated in [T81TT9] . The method allow to solve numerically nuclear- 
cascade equations for an arbitrary primary spectrum and for the most general 
form of hadron production cross sections. Here we do not focus upon de- 
tailed comparison with the results of other calculations, referring readers to 
the works [20lf2T1l22lf23ll2^1l25lf2^ and reviews [27p8lf29lf30] where one can find 
the relevant experimental data and comparison of calculations performed by 
other authors and a comprehensive list of references. 

The plan of this work is as follows. In sections 2 and 3 we touch briefly the pri- 
mary cosmic-ray spectra and modern high-energy hadronic models. In section 
4 we describe the method to solve the problem of meson transport through 
the atmosphere. In section 5 we calculate the hadron fluxes and compare them 
to old and new measurements having particular interest to test the calcula- 
tion scheme. In sections 6 and 7 we discuss the impact of hadronic models 
and the primary spectra on the sea-level muon flux at different zenith angles. 
In the short summary, we touch upon the uncertainties of current muon flux 
predictions resulted from the hadronic models and primary spectra. 



2 Primary cosmic ray spectra 

In our calculations we put primary emphasis on recent data on the primary 
cosmic ray (PCR) spectra and composition obtained with Advanced Thin Ion- 
ization Calorimeter, balloon-borne experiment ATIC-2 [8]. In order to compare 
the predictions with the high-energy measurements of the AM flux we extend 
the calculations to higher energies, up to 100 TeV, using also the PCR spec- 
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trum data of the GAMMA experiment [9]. The PCR model by Zatsepin and 
Sokolskaya [10], supported by the ATIC-2 data, was applied as the reliable 
instrument to extrapolate median-energy data to high-energy ones. Besides 
we use the parameterization of PCR spectra by Gaisser and Honda [TT] . 



The ATIC-2 experiment designed to measure primary cosmic ray spectra with 
the individual charge resolution from proton to iron in the wide energy interval 
50 GeV - 200 TeV enabled to obtain the data with a high statistical assurance. 
Proton and helium spectra measured in the ATIC-2 experiment have different 
slopes and differ from a simple power law (Fig. [I]). In this figure shown are also 
the energy spectra of protons and helium nuclei obtained in the balloon, satel- 
lite and ground based experiments 



The ATIC-2 data agree with the data of magnetic spectrometers (BESS [3T] . 
AMS [32], IMAX [33], CAPRICE [34j, MASS [35]) up to 300 GeV. In the 
energy region 1 < E < 10 TeV, the ATIC-2 data are consistent with the 
SOKOL [40J measurements and those of atmospheric Cherenkov light detector 
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Fig. 1. Primary proton and helium spectra combining balloon, satellite and 
ground-based measurements. The solid curves present the ZS model |10j . dashed 
- the GH parameterization [11] and dotted line - BV calculations 
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though the agreement is not so clear. 



The KASCADE experiment data for proton spectrum [41J in the energy range 
3 x 10 2 -10 6 GeV are deduced from the single-hadron (SH) spectrum assuming 
contributions from helium and heavy nuclei subtracted (green upward triangles 
in Fig. [1]). The later KASCADE measurements [12] for elemental groups of 
cosmic rays via the detection of extensive air showers (EAS) at energies above 
10 6 GeV (near and above the knee in the primary spectrum) are presented 
here by the primary proton and helium energy spectra (extracted directly from 
Figs. 14, 15 in Ref. [12] )■ These spectra were obtained on the base of EAS 
simulations with two hadronic interaction models, QGSJET01 (solid upward 
triangles) and SIBYLL 2.1 (solid downward triangles). 

Solid curves in Fig. [T] present the model suggested by Zatsepin and Sokolskaya 
(ZS) [10J that fits well the ATIC-2 experimental data and describes PCR spec- 
tra in the energy range 10-10 7 GeV. The dashed lines show the GH parame- 
terization [11]. In order to extend our calculation to higher energies, the PCR 
spectra measured in the GAMMA experiment was used. The energy spectra 
and elemental composition, obtained in the GAMMA experiment, cover the 
10 3 -10 5 TeV range (shaded bands) and agree with the corresponding extrap- 
olations of known balloon and satellite data at E > 10 3 TeV. In the present 
calculation employed was the version of the GAMMA spectra reconstructed in 
the framework of 1, 2D combined analysis with the SIBYLL 2.1 model (see [9] 
for details). 

The ZS model comprises contributions to the cosmic ray flux of three classes 
of astrophysical sources like supernova and nova shocks. In Fig. [1] also plotted 
are the PCR spectra (dotted lines) calculated by Berezhko & Volk (BV) [16] 
with use of nonlinear kinetic theory for diffusive acceleration of charged par- 
ticles in supernova remnants in which the magnetic field is amplified by the 
efficiently accelerating nuclear CR component. The BV model gives a natural 
explanation for the observed knee in the Galactic CR spectrum. 

The BV spectra and ZS ones are close in the energy interval 10-10 5 GeV 
inspite of the difference in shape, reflecting basic assumptions of the models. 
The difference in proton spectra becomes more appreciable at the energies 
10 5 - 10 6 while the difference in helium spectra is revealed above 10 6 GeV. 
Therefore one may expect that both ZS primary spectra and BV ones would 
lead to similar muon fluxes in wide energy range. 

At very high energies the BV spectra agree rather well with the KASCADE 
protons (QGSJET 01) and helium (SIBYLL 2.1) as well as protons and helium 
spectra derived from GAMMA experiment with usage of SIBYLL 2.1. Note 
also that some discrepancy exists between the KASCADE proton spectum 
and GAMMA one, both obtained with usage SIBYLL 2.1. 
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3 Hadronic interaction models 



To present day the direct measurements of inclusive cross sections for the nu- 
cleon and meson production in hadron-nucleus collisions have been extended 
up to energies of about 1 TeV only. Therefore cosmic-ray calculations require 
either developing a well-justified procedure for extrapolating cross sections 
measured at moderate energies or developing the models that would give re- 
liable predictions at high and ultrahigh energies. 

In our studies, we apply several hadronic interaction models, QGSJET 01 pQ, 
QGSJET 11-03 0, SIBYLL 2.1 @], NEXUS 3.97 0, and EPOS 0, that 
was successfully tested in the modern simulation programs CORSIKA |47j . 
AIRES [H] and CONEX [39] to analyze data on extensive air showers. Another 
hadronic model we use is the original one proposed by Kimel and Mokhov [TJ, 
for which we take updated parameters [T71I25II29] . Based on accelerator data, 
obtained at energies up to ~ 1 TeV, this model is compatible both in shape 
and magnitude with the current high energy hadronic models (Fig. [2]). The 
KM model was also applied in 3D Monte Carlo calculations of the atmospheric 
muon and neutrino fluxes [50J. 

Prominent features of /iA-interactions are the violation of Feynman scaling 
at high energy and the growth of total inelastic cross sections with energy 
(Fig. [3]). Figure [3] presents the hadron model predictions for the proton-air 
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Fig. 2. Spectra of the particle production in proton-air collisions at Eq = 10 TeV 
(left) and 100 TeV (right) according to the currently employed hadronic interaction 
models, KM, QGSJET, SIBYLL, NEXUS, and EPOS. 
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Fig. 3. Proton-air inelastic cross-section: the model predictions 
vs. experiments |51|52|53|54|55|56|57j . Solid line presents the fit 
<Jp\ = 290 — 8.7 log E + 1.14 log 2 E [EE]. 

inelastic cross section along with experiments [51 , 52.53, 54 55 56, 57J . Solid line 
here shows the fit to the data obtained with the prototype of the KASCADE 
hadron calorimeter [51~] . 

As illustrated in Fig. [2j the current models differ in the predicted amount 
of the scaling violation in the energy range 10-100 TeV. For example, in the 
central region (x < 0.1) QGSJET II predicts the growth of the inclusive 
cross section for reactions pA — > pX by a factor of 1.4 and for EPOS scaling 
violation is about 1.5. Apparently, this do not significantly affect the numerical 
results, because the contribution of small x to the fluxes of secondary hadrons 
is negligible due to falling energy spectrum of primary cosmic rays. Unlike 
the central region, the scaling violation in the fragmentation region (say, at 
x ~ 0.5 — 0.9) is most important, in particular in case of EPOS model: this one 
predicts a decrease of inclusive cross sections for pA — > pX up to 15%, while 
the EPOS pion and kaon yield is sizeably amplified (see Fig. [2] and Tabled]). 

To illustrate the difference of above hadron models the spectrum-weighted mo- 
ments was computed for proton-air interactions pA — > cX with the inclusive 
spectra F pc (E ,E c ) = (x/ o p n A ) da pc / dx: 




(1) 
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Table 1 



Spectrum -weighted moments Zp C calculated for 7 = 1.7. 



Model 


Eo, TeV 


z pp 


Zpn 




z pK+ + z pK~ 


QGSJET II 


0.1 


0.174 


0.088 


0.075 


0.0062 




1 


0.198 


0.094 


0.064 


0.0061 




10 


0.205 


0.090 


0.059 


0.0061 


EPOS 


0.1 


0.199 


0.112 


0.058 


0.0067 




1 


0.167 


0.084 


0.083 


0.0099 




10 


0.147 


0.073 


0.083 


0.0139 


SIBYLL 2.1 


0.1 


0.211 


0.059 


0.068 


0.0148 




1 


0.209 


0.045 


0.070 


0.0143 




10 


0.203 


0.043 


0.068 


0.0124 


KM 


0.1 


z pp 

0.178 


0.060 


Zp n + z pir~ 

0.044 0.0027 


z pK+ z pK- 

0.0051 0.0015 




1 


0.190 


0.060 


0.046 0.0028 


0.0052 0.0015 




10 


0.182 


0.052 


0.046 0.0029 


0.0052 0.0015 


where x = E c f E , c 


= P,n, 


TT ± ,K ± . 


As one can see 


from the Table 



demonstrate approximate scaling law of SIBYLL 2.1 and KM models, while 
in the case of QGSJET II and EPOS a sizeable violation of the scaling is 
found, particularly for the latter (see also [3|23] ). 



4 Transport of cosmic-ray mesons 



Here we concentrate on the meson component of the atmospheric hadron cas- 
cade, since the nucleon component was studied in detail [T7j. At sufficiently 
high energies, there is no need to consider 3D effects of the cascade, the electro- 
magnetic energy loss of hadrons, and the effect of the geomagnetic field. The 
collisions of cosmic-ray nuclei are treated with the semisuperposition model 
(SSP) [5B] (see also [55]) according to which the interaction of a nucleus (A, Z) 
with energy is similar to the averaged superposition of Z proton and A — Z 
neutrons interactions each with energy E = Ea/A. The SSP model predicts 
the same average values of additive observable quantities as the superposition 
model but the fluctuations. 

The meson component of the cascade decouples from the nucleon component if 
one neglects the small contribution of NN pair production processes in meson- 
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nucleus interactions. Under these assumptions, a set of integro-differential 
equations for charged-pion transport can be written as 



dir ± (E, h, •&) _ tt^M) m^7r ± (E,/i,^) 
dh K{E) pr 7r p(/i,^) 

+ E Gfe (E, M)+£ GkU (E, h, 0) + (2) 

i K 

OO 

1 r 1 da n ± n ±(Eo,E) ± 
+ TTTn / , in 7 F x TB 71 {E ,h,$)dE + 



K{E) J a? A (E) dE 



E 

oo 



K(E) J. a™ A {E) dE 



m i n 



where 7r ± (i?, /i, ^) is the flux of charged pions with energy E at the slant 
depth h (in gem -2 ) in the atmosphere along a direction of zenith angle 9, 
a™ A (E) is the inelastic-interaction cross section, X n (E) = [Nocr^^E)]" 1 is the 
pion interaction length (iVo = N A /A), p, m^, and are momentum, mass 
and lifetime of the pion, da(E , E)/dE is the cross section of charged pion 
production in ZiA-collisions, p(h, 9) is the air density profile, parameterized 
according to the Linsley's model of the US standard atmosphere (see [H] for 
details). The initial conditions for Eqs. (j2J) are ^(E, h = 0, 9) = 0. 

Sources of pions are strong interactions i+A = n +X (i = p, n, K ± , K°, K°) 
and the weak decays of kaons, K = K ± , K L , and K s . Hence two types of the 
pion source functions are 

OO 



E 



mm 



G%±{E,h,*)=B(K aK ) - ™* J ^F^(Eo,E)K(E ,h,$) (4) 



T K p(h,$) / p\ 



k 2tt 



j^raax 
K 13 



Here B(K 27T ) and B(K® 3 ) are the branching ratios for the decays K ± — * n^ir , 
Kg — > 7r + 7r - , and — > n^i^Ut. E , p , m K , and are the kaon energy, 
momentum, mass, and lifetime. and _F£ are the 7r-meson production 

spectra in the two- and three-particle kaon decays (see [60|61] ). The function 
Di(E Q , hjfl) in the right-hand side of (jHJ) and the function K(E ,h,'d) in (Hj) 
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are the fluxes of the parent particles i and K. The differential nucleon fluxes 
D P (E , h) and D n (E , h) were calculated using the formulas from the work [TT] . 
The processes of the pion regeneration and inelastic charge exchange 
counted by the last two terms in Eqs. (j2J). The integration limits in Eqs. (J2J), 
(J3J), and (j3J) are given in [T8] . 

At the first step we neglect a small contribution of kaons to the pion flux. This 
will be done at the second step after calculation of the kaon component. Let's 
denote by ^ the pion flux in the absence of the kaon source and introduce 
the combinations 

n ± (E, h, ■&) = tx + {E, h, 0)± Tt~(E, h, 0), (5) 
which obey the equations 



dIl ± (E, h, $) _ II ± (E, h, ■&) mJT^M) 



dh 



+ G% n (E } h } $) + 



(6) 



where x = E/E and 
G% K (E,h,#) = \G^ + (E,h,#) + G i % + (E,h,#) 



± 



Gg-{E,h,#) + G%-(E,h,#) 



E 



<a{E) 



da 7T + 7T +(E , £■) _^ da^+ w -(E Q , E) 



dE 



dE 



E Q =E/x 



Following [T7] we introduce the basic function of the method, ^-factor, which 
depends on the variables E, h, i? 



± 

Z±(E,M) = J $t(E,x 



IL^E/Xthtd) dx 
IL±(E,h,&) 



(7) 



and recast Eq. 



pr n p(h, i?) 



+ G^(£;,M) + 



The function /i, i?) serves as a generalization of the CR spectrum- weighted 
moments ^1/(7) (see e.g. [2311241162] ) to the case of non-power cosmic ray spec- 
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trum, violation of Feynman scaling, and energy- dependent inelastic cross sec- 
tions. 



Formally a solution of Eq. (jSJ) is: 



M) = / dt G% n (E,t,#) 



x exp 



I A 7r (,E) pr 7 ,p(z,d) 



(9) 



Then we solve via iterations the equations with unknown ^-factors. In zero- 
order approximation, we neglect the pion regeneration and charge exchange, 
Z±^(E,h) = 0. The differential energy spectrum of pions calculated under 
this assumption 



h r h / 

IT^ ) {E,h,ti) = I dt G% n (E, t, ■&) exp - J dz I — 



+ 



(10) 



enables one to obtain -2-factors in the first-order approximation. In the first- 
order and next approximations, the regeneration and charge exchange pro- 
cesses are already involved: 



H ± W(E,h,0) = Jdt G% ir (E,t,$) 



x exp 



dz ( l-Z% n \E,z,ti) + m n 
\ K{E) pT n p(z,$) 



In the (n+l)-iteration, the value of Z is calculated using the pion flux, ob- 
tained at the n-th step: 



± 

Z^ +1 \E,h,$)= J ^(E,x) 



n ±{ ri(E/x,h,ti) 



dx. 



(12) 



Now the process of iViV pair production in pion-nucleus collisions can be taken 
into account as a correction to the nucleon flux. In turn, the pion flux is refined 
taking into consideration this additional source. 



The evolution equations for charged and neutral kaons in the atmosphere have 
the form 
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dK(E,h,$)_ K(E,h,#) m K K(E,h,$) , /, tf) + 



dh ^k(E) pT K p(h,-d) 



'.Ml ^^^ K{W)iE °< (13) 

rrrnin 
KK 

where Gik{E, h, i?) are the source functions of kaons produced in NA- and 
7rA-interactions: 



pmio 

These equations are similar to Eqs.(E]), therefore we solve them using the same 
procedure. Finally we have 



KM(E,h,&) = J dt ^2 G iK (E, t, i?) x 



V,7T 



x exp 



d2 



+ 



pr K p(z, t?) 



(15) 



At last, we use the solutions, obtained for the kaon spectra, to take into 
account the contribution to the pion flux of the kaon source. 



5 Hadron fluxes 



The solution of the nucleon and meson transport equations discussed above 
enables us to compute the hadron fluxes with hadronic models under consider- 
ation and compare them to experiments. Fig. H] shows the differential intensity 
of nucleons in the atmosphere calculated for the GH primary spectrum along 
with the measurement data [63|I64|65| at depths of 20, 200, and 530 gem -2 . 
Evidently at small atmospheric depths there is no strong dependence of the 
flux on the hadronic interaction model but the difference accumulates as depth 
increases. The data shown in Fig. H] being obtained with rather large exper- 
imental errors do not afford a possibility to distinguish between this work 
predictions. The calculations performed with QGSJET II result in the maxi- 
mum nucleon flux whereas EPOS leads the minimum one. The rest of hadronic 
models spread between them. 

Experimental data on the neutron component at the sea level are more detailed 
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Fig. 4. Nucleon energy spectra at three atmospheric depths, computed with the GH 
primary spectrum and various hadronic models. The experimental data are from 
Refs. [63164165] . 



but rather contradictory. In Fig. [51 we compare the calculated energy spectra 
of neutrons at the sea level with the data taken from ref. [66] • The calculations 
with QGSJET01, KM and SIBYLL models result in best agreement with most 
recent data obtained from a hadron calorimeter prototype of the KASCADE 
facility [66]. One may hope that further experiments to study the high-energy 
nucleon flux will allow a more detailed test of hadronic models. 



Fig. [6] shows the calculated spectra of hadrons compared with the measure- 
ments [66 67 68 69 7rjl[Tip72|l73ll74|l75j . The hadron fluxes are computed at sev- 
eral depths of the atmosphere in the range 60-1030 gem -2 using the QGSJET 
II hadronic model, the ATIC-2 data, and the ZS cosmic ray spectrum in the 
region of higher energies (solid lines). Predictions for the GH parameterization 
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Fig. 5. Neutron energy spectra at sea level. The curves present this work calculations 
with usage of diverse hadronic models and the GH primary cosmic-ray spectrum. 
The experimental data are taken from Ref. |66j . 

combined with the QGSJET01, QGSJET II, SIBYLLL, NEXUS, and EPOS 
models are shown as shaded areas. Upper bounds of the bands correspond to 
the predictions of QGSJET II and lower ones - to those of EPOS, whereas 
the other models give the hadron intensities between the above predictions. 
Calculation results are scaled with the factors shown at the left side of the 
figure. 

At shallow depths there are only few experimental data sets available. For 
a depth of 60 gem -2 , we compare our results with the data of the Moscow 
group [67] obtained onboard balloons with X-ray emulsion chambers. For a 
depth of 260 gem -2 we use the data [68] (shown as rectangle area) obtained 
on the airplane with help of emulsion chamber also. In both the cases, our 



13 



> 
o 



05 



10 
1 



lO" 1 ! 



10 
10 



-2 



■3 ^ 



10" 4 
10 

T3 



10 
10 
10" 
10" 
10 



10 



10 11 
lO" 12 



I I I I 1 1 1 1 — I — I I I I 1 1 r~ 



-rrr| 1 | 



I 

h, g/cm 
60 



260 




a Pamir (2005) 

t Kuzmichev et al. (1983) 

* Thien-Shan (1983, 89) 

Erofeeva et al. (1968) 

& Babayan et al. (1967) 



_i_ 



• EAS-TOP (2002) 

A KASCADE (1995) 

□ Shmeleva et al. (1983) 

o Ashton etal. (1975) 

■ Brooke et al. (1964) 

I I I I I I I i i i 



10 2 



10 J 



10 4 



10 5 10 6 
E, GeV 



Fig. 6. All hadron energy spectra in the atmosphere. Solid lines present this work 
calculation with the QGSJET II hadronic model and ATIC-2 primary cosmic-ray 
spectrum. The bands show the spread of our predictions due to the usage of diverse 
hadronic models for the GH primary spectrum. The experimental data are from 
Refs. |66l67l68l69l7()l7ll72l73l74l75] . 

results well agree with the data within the experimental errors. 

For atmospheric depths of 600-700 gem -2 , we compare our predictions with 
measurements performed at the Thien-Shan and Pamir. High-energy data ob- 
tained with the carbon emulsion chambers in the Pamir experiment [69] agree 
well with our calculations, while the data from the ionization calorimeter of 
the Thien-Shan experiment [70] give the hadron intensity at E > 50 TeV 
about twice as high as our prediction with QGSJET II model. For the atmo- 
spheric depth 820 gem -2 , our calculations agree rather well with the EASTOP 
data [73j. For 1030 gem -2 our calculations at energies above TeV are in some 
discrepancy with the data of the KASCADE experiment |66j as well as with 
those of early measurements (taken from [66J). Also in Fig. [TJwe compare in 
detail our predictions for each of considered hadronic models with the recent 
data of the Pamir and EAS-TOP experiments. 



In Fig. El we plot the calculated ratio of the pion flux to nucleon one (n/N 
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Fig. 7. Scaled data of the Pamir and EAS-TOP experiments. Lines present the 
calculations for the GH primary spectra using every hadronic model under study. 

ratio) at sea level together with experimental results and other calculations 
taken from [66J. As one can see, the primary spectra and mass composition 
weakly affect the ratio unlike hadronic models. This allows in some way to 
test the cross sections: the models KM, NEXUS, QGSJET 01, and QGSJET 
II yield comparatively similar results, whereas the models SIBYLL and EPOS 
lead to the sharp difference in tt/N ratio. The left panel of Fig. [9] shows the 
charged K/tt ratio of the fluxes calculated for the atmospheric depth 200 
gem -2 . It is seen that the NEXUS and EPOS predict maximal K/tt ratios 
~ 0.2-0.25 at E = 100 TeV, whereas the KM and QGSJET II give minimal 
values, ~ 0.1-0.15. The right panel of Fig. [9] shows the same ratio unfolded 
along the atmospheric depth. 



6 Muon fluxes 

The muons are abundantly generated in decays of the mesons that are pro- 
duced through hadron-nuclei collisions in the atmosphere. Apart from the 
dominant source of atmospheric muons, tt^ and K^ decays, we take into 
consideration three-particle semileptonic decays, K^ 3 , K%, and small con- 
tributions originated from decay chains K — > tt — > // (Kg — > tt + + 7r~, 
K^ — > ^ + 7T°, K\ — > n ± + £ T + Pifa), £ = e, //). At very high energies, the 
bulk of muons is expected to arise from decays of the charmed particles, D 
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Fig. 8. The ratio 7r/(p+n) calculated at sea level with the set of hadronic models. 
The ATIC-2 primary spectrum and GH parameterization are used. Experimental 
data (closed symbols) and calculations of other works (open symbols) are taken 
from Ref. |66j . 
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Fig. 9. Charged i^/7r-ratio, depending on the energy, calculated for the atmospheric 
depth of 200 gcm~ 2 (left plot), and that at 10 TeV unfolded along the atmospheric 
depth (right plot). GH primary spectrum is used. 
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and A c . 



The results of this work calculations of the surface muon flux near the vertical 
are presented in Figs. [10] and [11] along with data of the old and recent muon 
experiments: the direct measurements of the CAPRICE [12], BESS-TeV [13], 
L3+Cosmic [14J, Cosmo-ALEPH (see Ref. [15]), MASS [76], and L3 (taken 
from [21]) as well as data (converted to the surface) of the underground ex- 
periments MSU [77], MACRO [78], LVD [79] , Frejus [80], Baksan [gT], and 
Artyomovsk [82J. 

The light shaded area in these figures shows the muon flux computed with 
the KM hadronic model spectrum taking into consideration the uncertainties 
in the ATIC-2 spectrum data. The dark shaded area on the right in Fig. [TT1 
shows a high-energy part of the muon flux calculated with the KM model for 
the GAMMA primary spectrum. The curves represent muon flux predictions 
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Fig. 10. Energy spectrum of muons at ground level near the vertical. The curves 
show this work calculations for the ATIC-2 primary cosmic-ray spectrum and dif- 
ferent hadronic models. The band presents a result for the KM model taking into 
consideration statistical errors of the ATIC-2 data. The experimental data are 
from [12113114115124176177178179180181182] . 
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Fig. 11. High-energy plot of the ground level muon spectrum. The solid curves 
and bands show this work calculations for the KM model with usage of ATIC-2, 
GAMMA, and ZS primary spectra. The reults for SIBYLL (dash-dotted line), 
QGSJET II (dashed), and EPOS (bold dotted) are obtained using the ZS primary 
spectrum. The dotted lines present sum of the conventional and prompt muon flux 
(see text for details). 

we have made basing on the ATIC-2 primary spectra and using the variety 
of hadronic models: solid, thin, dashed, and bold-dotted lines stand for the 
calculation with usage of the KM, SIBYLL, QGSJET II, and EPOS models 
respectively. Note that at energies beyond the scope of the ATIC-2 experiment, 
the ZS primary spectrum model was used. The latter fits well the ATIC-2 data 
and can be considered as a reliable bridge from the TeV to PeV energy range. 

The muon flux predictions obtained with usage of the KM, SIBYLL, and 
EPOS models agree well with recent experimental data for energies below 10 
TeV. The range of the muon flux predictions indicates the minimal uncertainty 
(~ 10%) due to hadronic models . This uncertainty is comparable with that 
in the primary cosmic ray spectrum measured in the ATIC-2 experiment. 
However, in a case of the QGSJET II model the muon flux is 30% lower if 
compared both with the experimental data and calculations mentioned above. 
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The predictions of the AM flux with usage of the GH primary spectrum were 
obtained for the KM hadronic model in Ref. [18] and for the QGSJET 01, 
SIBYLL, NEXUS, and QGSJET II models in Ref. [S3]. 

The prompt muon component of the flux is shown in Fig{TT] (dotted lines) 
where numbers above the curves indicate the calculations based on the charm 
production models: 1 - quark-gluon string |84.85j , 2 - recombination quark- 
parton [24"|85] . and 3 - Volkova et al. [86J. These components were added as 
corrections to the conventional flux (solid line) in order to demonstrate the 
extent of the uncertainty of the calculated muon flux at PeV scale. 



7 Zenith-angle distribution of the muon flux at sea level 

The zenith-angle distribution of the atmospheric muon is of a special interest 
because it is sensitive to details of the meson production and decays [20]. In 
particular, the K/n flux ratio, being dependent both on the hadron model 
and primary cosmic ray composition, has effect on the muon flux at inclined 
directions. The consistent description of the flux at large zenith angles and at 
the vertical direction would confirm the adequacy of the hadronic model and 
the validity of the primary cosmic ray spectrum data. 

A comparison of the measurements in wide zenith angle range with the muon 
energy spectra calculated for the ATIC-2 primary cosmic ray spectra and KM 
model is shown in Figs. [T2"] - fT6l The calculations agree rather well with old ex- 
perimental data obtained with magnetic spectrometers, AMH [87], BMS |88j . 
Kiel-DESY [S3], MUTRON [21], DEIS [HI] or X-ray emulsion chambers of 
MSU [T7t9"2"] . as well as with comparatively recent results obtained by Okayama 
telescope [93], those of the Karlsruhe experiment [91] and the L3+Cosmic |14j . 

Figure [H] shows the reliable measurements with the magnetic spectrometers 
DEIS and MUTRON for zenith angles ranging from 78° to 90° and our predic- 
tions based on the KM model(shaded areas) and QGSJET II (hatched area). 
For clearness, the QGSJET II predictions are shown only for the two angle 
ranges 78°-80° and 89°-90°. Unlike the vertical direction, the muon data for 
large zenith angles are described by QGSJET II rather well. A difference be- 
tween the QGSJET II and KM models reduces now to 15-20%. Muon fluxes at 
angles of 45°, 72°, and 89° were measured by the MSU group [77|92] using with 
X-ray emulsion chambers (see Fig.[T4" l[l~5"|) . The calculations done using the ZS 
primary cosmic-ray spectrum and KM model agree well with the MSU data 
including horizontal directions. It is evident from Fig. [15] that the difference 
between the QGSJET II and KM models decreases approximately from 30% 
to 15% as zenith angle increases. To compare the predictions with the recent 
L3+Cosmic experiment data for the zenith-angle range 0° < 9 < 58° (Fig. fTB]) . 
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p, GeV/c 

Fig. 12. The muon flux calculated for wide range of zenith angles in comparison with 
the data of AMH [87] magnetic spectrometer. Calculations are performed for the 
ATIC-2 primary cosmic ray spectra and KM model. The bands show the calculated 
flux spread due to uncertainties in zenith angles. 

we apply the atmospheric density profile measured near the site of L3+Cosmic 
experiment. Apparently our calculation agrees with the L3+Cosmic data at 
Pfi ~ 50 GeV/c. Solid curves correspond to the average zenith angles for a 
given angle bin. Also for comparison, the low-energy data of the MASS exper- 
iment [76] are added to the figure. 

The effect of Kj it ratio on the zenith-angle distribution of high-energy muons 
is rather transparent: each new source of the muon flux leads to decrease of the 
muon flux ratio, D tl (9)/D tl (0) [201195] . The more intensive is source (i.e. the 
larger K/ir ratio in this case), the lesser is the angle enhancement of the muon 
flux. However we can observe the influence of K/tc ratio (that is the proper 
feature of the hadron model) only in the restricted energy region, because the 
ratio D fJj (9)/D /M (0) weakly depends on the K/tt ratio for E 3> where 
€k(@) is the critical energy for kaons (e.g., e#±(0) ~ 850GeV). The reason 
is that the meson decay probability decreasing with the energy, Wd(E, 9) oc 
e^ K {9)/E, leads to D^/D^O) oc e^ K (9)/e n , K (0), or D^/D^O) ~ 7.5 
near the maximum of the muon zenith-angle distribution (see Fig. [T7|) . 

Fig. [T71 shows the muon flux ratio D^(9)/D^(0) as a function of cos#, that is a 
coefficient of zenith-angle enhancement of the muon flux, calculated with the 
EPOS model (solid lines) and QSJET II one (dashed) in the range 1-50 TeV. 
A choice of these hadronic models is motivated by an observation that the 
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Fig. 13. Comparison of predicted muon fluxes for inclined directions with the data 
of Brookhaven [88], Kiel-DESY [89], MSU [77192] . Okayama [93] and Karlsruhe 
experiments. Same notations as in Fig. fl2l 



EPOS and QGSJET II lead to the outermost predictions of the hadron fluxes 
and K/ir ratio (see Fig. [7J and [9]). This gives a possibility to examine indirectly 
the effect of the kaon source on the muon angular flux: the EPOS prediction 
of the K/ir ratio (~ 0.2 at E = 10 TeV) results in lower D fl (6)/D fl (0) ratio in 
comparison with that of QGSJET II model for which the K/ir ratio is close 
to 0.1 at the same energy (compare also with predictions by Volkova [H2] for 
K/tt = 0.1, closed circles in Fig. |T7|) . 

A comparison of the calculated ratio with measurements in a wide en- 

ergy range makes a possibility to study hadron-nucleus interactions indirectly. 
Present computation (solid line) of the muon charge ratio at sea level is 

plotted in Fig.[T8]along with experimental data from Refs. [96l97i98 U99|ll00|ll01|ll02l 
This computation has been performed with the KM hadronic model and the 
ATIC-2 primary spectrum. For muon energies below the charged pion critical 
energy (~ 115 GeV), the pions are the mostly dominant source of muons. 
For these energies, the calculated ratio, ~ 1.3, agrees with the recent data of 
BESS-TeV, CosmoALEPH, and L3+Cosmic experiments. As energy increases 
the kaon source of muons intensifies leading to a maximum of the muon ra- 
tio, ~ 1.4, at energy close to 10 TeV. This value agrees with the most recent 
results obtained with the MINOS far detector [101] . 
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p, GeV/c 

Fig. 14. Comparison of near-horizontal muon fluxes at sea level with the data of 
MUTRON [90] and DEIS [91J. The calculations are made for the ATIC-2 primary 
spectrum with usage of KM (shaded areas) and QGSJET II (hatched areas). 

8 Summary 

In this work, we have calculated the atmospheric high-energy fluxes of nu- 
cleons, mesons and muons implementing the hadronic models QGSJET [Tp] . 
SIBYLL @], NEXUS [5J, EPOS [6], and the Kimel and Mokhov one [7]. 
These calculations basing on the method [PTfPcS] have been performed using di- 
rectly data for primary spectra and composition of the ATIC-2 experiment [8], 
GAMMA one [9] and Zatsepin and Sokolskaya model [10]. Alternatively we 
used well known parameterization by Gaisser and Honda [11]. The method ap- 
pears rather efficient tool to study the transport of cosmic ray particles through 
the atmosphere in a generalized case of the non-power law primary spectrum 
and non-scaling behavior of hadronic cross sections. Various hadronic models 
can be rather easily embedded into the developed code. 

The present calculations demonstrate the consistency of the new primary 
comic ray spectra measurements with recent experimental data on atmospheric 
hadron and muon fluxes at different levels of atmosphere for a wide range of 
energies and zenith angles. The comparison of the calculated and measured 
fluxes gives the possibility to evaluate the uncertainties originated from the 
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Fig. 15. Comparison of the oblique muon fluxes at sea level with the data of MSU 
experiments |77|92j (see also Fig. [T3|) . The calculations are made for the ATIC-2 
primary spectrum using KM model (solid lines) and QGSJET II (dashed). 
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Fig. 16. Comparison of the predictions for near-vertical muon fluxes with the data 
of the L3+Cosmic experiment. The computations are made for the ATIC-2 primary 
spectrum with usage of KM. The MASS data are added to view the low-energy 
junction. 
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Fig. 17. Zenith-angle enhancement of the muon flux, D^(E,9)/D^(E,0°), calcu- 
lated with the EPOS (solid lines) and QSJET II (dashed), both for the ATIC-2 
primary spectrum. Symbols present the calculations of [22] (squares) and [95] (cir- 
cles). Numbers above curves are the muon energies in TeV. 

errors of the primary spectrum data. The high accuracy of the ATIC-2 data 
resulted in the muon flux calculation uncertainty, comparable with rather high 
precision of the last decade muon flux measurements. Also this allowed us to 
study the effect of the non-power law structure of the ATIC-2 primary spec- 
trum on secondary fluxes. The effect on the muon spectrum shape of deviations 
from a power law in Zatsepin and Sokolskaya primary spectrum appears at 
energy above 10 TeV. 

The vertical muon flux predictions, obtained with usage of the ATIC-2 primary 
spectrum and KM, SIBYLL, EPOS hadronic models, are in accordance with 
recent experimental data in the wide energy range. These models give similar 
muon fluxes differing by ~ 10%, that is comparable to the uncertainty in 
the ATIC-2 data. The muon spectrum predicted with the QGSJET II model 
is ~ 30% lower in the energy range 10 2 — 10 5 GeV if compared both with 
the experimental data and calculations with usage of the rest models. For 
near horizontal directions the uncertainties in the muon flux calculations are 
found to be less in comparison with that for the vertical. In particular, the 
discrepancy between the QGSJET II and KM models is reduced to something 
like 15%. Thus we may state: the uncertainties due to hadronic cross sections 
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Fig. 18. Muon charge ratio at sea level calculated with KM model for 
the ATIC-2 primary spectra. The data of experiments are taken from 
Refs. |13I14I96I97I98I99I100I101TT02] . 



are about 30%, while those from variations in the primary spectrum are at a 
level of 10%. 



And at last, the prompt muon component originated from the decay of the 
charmed hadrons is still not be reliably extracted from the experiments, and 
it seems to be the source of hardly estimated uncertainty that superimposes 
on the uncertainty factor due to poor knowledge of the primary cosmic ray 
flux around the knee. 
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